rm(list = ls(all = TRUE))
setwd("C:\\Users\\lenovo\\Desktop\\R")
#13.1 二分类结果的评价指标
#将某一检测指标的测定值与金标准（gold standard）相比较
#金标准为特异性诊断
#如果没有金标准，也可以用医学专家共同制定的公认的综合诊断

#13.1.1 灵敏度与特异度
#灵敏度（sensitivity），真阳性
#患者检测结果为阳性的概率
#特异度（specificity）。真阴性
#非患者检测结果为阴性的概率

#13.1.2 预测值
#阳性预测值（positive predictive value）
#在检测结果为阳性情况下受试者患病的比例
#阴性预测值（negative predictive value）
#在检测结果为阴性情况下受试者不患病的比例

#13.1.3 图解灵敏度、特异度、阳性预测值和阴性预测值
table1 <- as.table(cbind(c(80,20),c(10,90)))
dimnames(table1) <- list("检测结果" = c("阳性","阴性"),
                         "金标准" = c("有病","无病"))  
table1
mosaicplot(t(table1), col = c("red","white"), main = "")

table2 <- as.table(cbind(c(80,20),c(190,1710)))
dimnames(table2) <- list("检测结果" = c("阳性","阴性"),
                         "金标准" = c("有病","无病"))  
table2
mosaicplot(t(table2), col = c("red","white"), main = "")

#患病率很低的筛选检测中很少用到阴性预测值

#13.1.4 诊断试验的综合评价指标
#1 正确率
#正确区分患者、非患者的能力
#但在很大程度上依赖受试人群的患病率

#2 约登指数（Youden Index）
#灵敏度+特异度-1
#在-1到1之间，
#值越接近1，诊断试验真实性越好
#值为负数，没有临床应用价值

#3 似然比
#阳性似然比（positive likelihood ratio）
#真阳性率/假阳性率 = 灵敏度/（1-特异度）
#阴性似然比（nagative likelihood ratio）
#假阴性率/真阴性率 = （1-灵敏度）/特异度

#似然比为1，说明在患者和非患者中得到该检验结果的概率相同，没有诊断价值


#13.2 ROC及曲线下面积
#ROC: Receiver Operation Curve
#AUC: Area Under Curve
#真阳性率（灵敏度）为纵坐标，假阳性率（误诊率）为横坐标

#13.2.1 单个ROC评价
#install.packages("pROC")
library(pROC)
data(aSAH)
str(aSAH)

roc1 <- roc(outcome ~ s100b, data = aSAH)
attributes(roc1)
roc1$auc

#求约登指数
roc.result <- data.frame(threshold = roc1$thresholds,
                         sensitivity = roc1$sensitivities,
                         specificity = roc1$specificities)
roc.result$youden <- roc.result$sensitivity + roc.result$specificity - 1
head(roc.result)

#找出约登指数最大的一行
which.max(roc.result$youden)
roc.result[18, ]

#可以一步完成
library(pROC)
#best 最佳截断点
#best.method 默认为"youden"
#也可以设为"closest.topleft"，选取ROC曲线最左上角
coords(roc1, "best", transpose = FALSE)
coords(roc1, "best", best.method = "closest.topleft", transpose = FALSE)

#绘图
plot(1-roc1$specificities, roc1$sensitivities, type = "l", lwd = 2)

library(pROC)
plot.roc(roc1, print.auc = TRUE, auc.polygon = TRUE,
         grid = c(0.1,0.2), grid.col = c("green","red"),
         auc.polygon.col = "lightblue", print.thres = TRUE)

#AUC的置信区间——DeLong法
ci.auc(roc1)

#13.2.2 两个ROC的比较
roc1 <- roc(aSAH$outcome, aSAH$s100b)
roc2 <- roc(aSAH$outcome, aSAH$ndka)

#DeLong非参数方法，Venkatraman回归模型法，bootstrap重抽样法
#默认为DeLong法
#默认为两组相关检测结果AUC的比较
#独立：参数paired = FALSE
roc.test(roc1,roc2)

#绘图
plot(roc1)
lines(roc2, col = "red")
test <- roc.test(roc1, roc2)
text(0.5,0.5, labels = paste("p-value = ",round(test$p.value, 3)))
legend("bottomright",
       legend = c("S100b", "NDKA"),
       col = c("1","red"), lwd = 2)

#13.2.3 Logistic回归的ROC曲线
#建立Logistic回归模型
fit <- glm(case ~ induced + spontaneous, family = binomial, data = infert)
library(epiDisplay)
logistic.display(fit)
#ROC绘制
lroc(fit, line.col = "red", lwd = 2)


#13.3 联合试验
#用两种或两种以上的诊断试验对疾病作出判断的方法称为联合实验
#平行试验（并联）：有阳则阳，全阴才阴：灵敏度高、特异度低
#系列试验（串联）：全阳才阳，有阴则因：灵敏度低、特异度高

library(pROC)
cut.point1 <- coords(roc1, "best", transpose = FALSE)$threshold
cut.point1

cut.point2 <- coords(roc2, "best", transpose = FALSE)$threshold
cut.point2

#平行试验
parallel <- ifelse(aSAH$s100b > cut.point1 | aSAH$ndka > cut.point2,
                   "Positive", "Negative")
table(parallel, aSAH$outcome)

#系列试验
serial <- ifelse(aSAH$s100b > cut.point1 & aSAH$ndka > cut.point2,
                   "Positive", "Negative")
table(serial, aSAH$outcome)

#降低漏诊率，用联合试验
#诊断慢性病、筛查疾病，用系列试验